1 Introduction

The proportional multi-state multiple-cohort life table model (PMSLT) is a population level model (macro) approach to simulate health (and economic) implications of changes in exposure to health risk factors (e.g. physical inactivity, air pollution and diet). The PMSLT has been widely used to simulate outcomes for population level interventions for the reduction of chronic diseases.

The model was developed by Jan Barendregt and colleagues and has been widely used in Australia and New Zealand (T. Vos et al. 2010; Blakely et al. 2015).

The basic infrastructure of the model consist of three components: (1) Effect size for the intervention of interest (e.g. intervention to urban design that modifies population levels of physical activity); (2) Calculation of the potential impact fraction (PIF) to derive the change in occurrence of disease (incidence rate/case fatality rate) attributable to a change in the distribution of the risk factor (e.g. physical activity); and (3) Use of the PMSLT to simulate health (and economic) outcomes attributable to a change in the distribution of health risk factor/s in the population of interest. Figure 1 summaries the basic infrastructure of the model. ITHIM is included in Figure 1 to show that both approaches share in common steps one and two and differ in the mechanisms of calculation of change in health burden.

HALYs, QALYs and DALYs

In this model we use the term health-adjusted life year (HALY). As summary measure of population health it measures both quantity and quality of life, where one HALY represent the equivalent of one year in full health (which could be two years with a quality of life of 0.5, for example). Specific types of HALY are the quality-adjusted life year (QALY) and the disability-adjusted life year (DALY). The QALY derives from economics and was first used in the 1960s as a measure of health gain (Gold, Stevenson, and Fryback 2002). The disability-adjusted life-year (DALY) was developed for use in burden of disease studies as a measure of health loss due to disease (Gold, Stevenson, and Fryback 2002). Our calculated HALYs are neither QALYs not DALYs, but something in between. They are similar to QALYs in that they represent health gains. However, the main difference is in the calculation of the health-related quality of life component. QALYs use measures of utility weights that traditionally represent individual experiences of health, whereas our estimated HALYs use disability weights linked to specific diseases, which were developed for the Global Burden of Disease study (Gold, Stevenson, and Fryback 2002). As discussed in past research (L. Cobiac, Vos, and Barendregt 2009; Roux, Pratt, and Tengs 2008) the main advantage of using disability weights over utility weights is that disability weights refer to specific diseases rather than health states (which are difficult to link to risk factors-e.g. physical inactivity). We opted to use the more general terms HALYs given that the use of the DALYs terminology may lead to think that our calculations are similar to those in burden of diseases studies (Murray et al. 2012). In our study, our model does not explicitly separate years of life lost (YLL) and years lived with disability (YLD) components, but instead calculates the total number of life years lived, adjusted for the average health-related quality of life in those years (by age and sex). In burden of disease studies, DALYs are defined as the sum Years of Life Lost (YLL) and Years Lived with Disability (YLD).

Basic ITHIMR infrastructure

Basic ITHIMR infrastructure

1.1 Contribution to ITHIMR

The PMSLT similar to ITHIM is a comparative risk assessment approach (Briggs, Scarborough, and Smith 2016) that consist of calculating the change in the health burden for a population of interest from a change in exposure to health risks factors (e.g. physical inactivity, air pollution and road trauma).As depicted in Figure 1, both methods need estimates of the potential impact fraction (PIF), which indicates the proportion of the disease burden attributable to a risk factor of interest (e.g. physical inactivity) (Barendregt and Veerman 2010). A step further back, is the development of scenarios that bring about change in the distribution of the risk factor of interest. For now, we only focus on calculations from the PIF onward, and provide a hypothetical example of change in the population distribution of physical activity. Incorporation of additional health risk factor (air pollution, road trauma, NO2 and noise) will be discussed in the relevant code sections.

1.1.1 Difference between ITHIM and PMSLT

  • Time component The PMSLT follows a population of interest over time. For example, as set up here, we simulate sex and age (5 years starting at 20) cohorts over time until they die or reach 100 years of age. This implies that we can include trends for diseases, time lags for change in exposure to risk factors and change in health and demographic changes (e.g. population growth). In addition, we can estimate yearly changes in the burden of diseases over the life course or for a specified number of years. The ITHIM approach is a snapshot of change in burden for one year.

  • Interaction between multiple diseases The PMSLT accounts for the interaction between multiple diseases, with proportions of the population being able to be in more than one health state (Briggs, Scarborough, and Smith 2016). This avoids overestimation of outcomes as a result of summing health outcomes attributable to each disease individually as done in ITHIM. It is important to note that the PMSLT assumes that diseases are independent of each other. That is to say, developing a disease is unrelated to a concurrent diagnoses of another disease).

  • Mortality rate The PMSLT calculations for changes in life years (and health-adjusted life years) and mortality outcomes is based on observed mortality rates for the population of interest. In the ITHIM model, if burden of disease estimates from the Global Burden of Disease (GBD) study are used, then, the mortality component is based on the highest attained life expectancy observed in the world.

  • Impact of disability in increased life expectancy In GBD studies, YLLs are not adjusted for disability; hence, their use in estimating intervention effects results in over-estimation, which the PMSLT approach avoids. Another way of seeing this is that estimated changes in morbidity using the ITHIM do not allow for how implicit increases in life expectancy impact on morbidity. While the changes in deaths and prevalence using the PMSLT are in some ways more accurate than those from the ITHIM approach it should be noted that that the average age of death and incident disease will change and thus the disease burden will be on average be shifted later in life (which is a realistic approach).

2 R development

The model is set up as a long script to perform the required mathematical calculations. Where possible, we wrote functions and loops to avoid repetition. We set up the model with data for Greater London. Figure 2 depicts the PMSLT model framework, which was followed in the code development.

Proportional multi-state life-table simplified framework. The simplied PMST shows the interaction between the life table, disease life table and potential impact fraction (PIF). The PIF calculations by age and sex group are the same as those generated for ITHIM. The PIF (or 1-PIF) modifies incidence of disease, which changes prevalence and mortality (disease specific life table). Changes in prevalence and mortality rates from the disease specific life tables feed into the life table by changing all-cause mortality, which in turn changes life years. Change in prevalence of diseases changes total years lived with disability, which in turn modifies health-adjusted life years

Proportional multi-state life-table simplified framework. The simplied PMST shows the interaction between the life table, disease life table and potential impact fraction (PIF). The PIF calculations by age and sex group are the same as those generated for ITHIM. The PIF (or 1-PIF) modifies incidence of disease, which changes prevalence and mortality (disease specific life table). Changes in prevalence and mortality rates from the disease specific life tables feed into the life table by changing all-cause mortality, which in turn changes life years. Change in prevalence of diseases changes total years lived with disability, which in turn modifies health-adjusted life years

In what follows, first, we specify input parameters. Second, we present the code with explaining notes. Third, we present examples of outcomes and lastly we comment on topics related to implementation. Here we only included the physical activity health pathway. In the comments section, we discuss the implementation of exposure to air pollution and road trauma. Note that in the presentation of input parameters, those needed to calculate PIFs are excluded, as these are common to the ITHIM, expect if trends are included (refer to comments section).

2.1 Inputs

We specify data requirements for the life table and disease life tables (Figure 2) and potential sources.

2.1.1 Life table

Inputs of the life table are: population numbers by sex (per 1-year or age grouping of interest), mortality rates or probability of all cause mortality by single age group and sex and total prevalent years lived with disability rate per single year by sex. Disease specific disability weights are presented as inputs here as these adjust the total years lived with disability, hence, the health-adjusted life years.

2.1.1.1 Population numbers

These data will be provided by the synthetic population or derived from other data. In the code presented here, we simulate 5-year age and sex cohorts. Data for population may be in five-year age groups or one-year. For the example for Greater London, we derive 5-year age groups from GBD IHME data, however, we also provide a script if five-year age groups are to be derived from one-year age groups data. I left potential data sources below as a reference.

Data source: (1) National census; (2) Worldwide population and mortality data: http://www.mortality.org/ (mostly high income countries; and (3) Calculate from the Global Burden of Disease by the Institute of Health Metrics and Evaluation (GBD IHME) data (rates and numbers available from (http://ghdx.healthdata.org/gbd-results-tool).

2.1.1.2 Mortality rates

Mortality rates are needed per single year and sex. These data are available from GBD IHME, however, in age groups (1-4, 5-9, etc). We provide an interpolation script to derive in between ages rates (cubic spline).

Note that we need data for population numbers and all cause mortality rates for: (1) PMSLT and (2) Dismod II collection (more in Dismod II section). Population data from the synthetic population is used for the PMSLT (if available). For Dismod II, population and mortality data should be from the same source (GBD IHME).

2.1.1.3 Total years lived with disability rates per single year and sex.

These data is available from the GBD (http://ghdx.healthdata.org/gbd-results-tool) per 5-year age groups. We can use interpolation to derive between ages rates (script provided).

2.1.1.4 Disability weights (quality of life weights)

Disability weights (DW) are derived from disease specific years lived with disability (YLD) and disease specific prevalence by age group (5 years) and sex. Data for YLDs prevalence is obtained from the online GBD IHME data tool (http://ghdx.healthdata.org/gbd-results-tool). An age and sex specific-correction was introduced to counteract the effects of accumulating comorbid illnesses in the older age groups (Equation 1).

\[\begin{equation} \label{DW adjusted for total YLDs-YLDdPd1YLDt} (YLDd/Pd)/(1-YLDt) = DW adjusted for total YLDs \end{equation}\]

Where YLDd is the YLD mean number per age and sex for a given disease, Pd is the prevalence (as reported in GBD ) for a given disease by age and sex and YLDt is total YLD rate per age and sex.

2.1.2 Disease life tables

2.1.2.1 Incidence and case fatality

For each of the modeled diseases the PMSLT needs incidence and case fatality rates per sex and one-year intervals. Data from the GBD IHME studies with Dismod II (free at https://www.epigear.com/index_files/dismod_ii.html) is used to derive internally consistent data and generate missing data. For example, the GBD studies provide data for incidence, prevalence and disease mortality, however, not case fatality. Other national level sources may also be explored/used, and compare with estimates produce from GBD data and Dismod II.

Dismod II inputs are: (1) population numbers and mortality rates and (2) disease specific inputs.

Population and mortality

Within Dismod II, each setting (e.g. country) has a collection that consists of population numbers (preferably the same as used in GBD IHME studies, due to the mortality envelop) and all- cause mortality rates (numbers and calculate rates). The GBD provides 5-year age groups that are acceptable input parameters for Dismod II.

Disease inputs by age group and sex

Each setting collection has a given number of diseases. Dismod II works with at least three of: case fatality, prevalence, incidence, mortality (disease), case fatality, remission, duration and the relative risk for mortality. So far, we have been assuming that remission is zero for chronic diseases, that is to say, when people become diseased, they do not recover. Special care should be taken with this assumption, as the GBD data assumes remission for some diseases, for example cancers, where after 10 years cases recover, except for long term sequelae. Since GBD now provides prevalence, incidence and mortality, it may be best to use all three as Dismod II input parameters to compare the effect of the remission assumption by the GBD for some diseases.

2.2 Code

Following the structure of Figure 2, we developed functions to perform sex and age cohorts calculations for the life table, disease life tables and potential impact fractions: run_life_table, run_disease and and run_pif. We also generated two functions for outputs: plot_outputs and gen_aggregate. The function plot_outputs creates age-group and sex linear plots for specified outcomes (e.g. health-adjusted life years, incidence of diabetes) and gen_aggregate adds up each cohort results. Functions were then used in a code script. In what follows, we explain each step in the development of the script. Here we also include code chunks, however, we also kept them separately in the MSLT folder, in the code file.

In what follows, we start explain the script step by step.

2.2.1 Set up

We start by cleaning the global environment (1) to keep track of our works and ensure that the code is generating our desired outcomes. Then, we set up an option to avoid the use of scientific notation (2) and lastly we load the functions (3). The code chunks are shown in the rmarkdown output.

  1. Clean Global Environment
rm (list = ls())
  1. Avoid scientific notation
options(scipen=999)
  1. Load functions
source("code/functions.R")

2.2.2 Inputs

Table 1 describes data needs for the PMSLT, here we expand on the data needs and mechanisms (Figure 3) to use the PMSLT approach in ITHIMR (Figure 1).

Initial case studies for the ITHIMR are: London, Sao Pablo, Delhi, Accra, Los Angeles and Edinburgh. Here, we will start with Greater London given the availability of disease epidemiology data from the GBD IHME study. For the rest of the case study cities data is available at the country level, hence, a scaling method is needed to reflect the local burden of disease.

Proportional multi-state life-table model. Three sections are presented in Figure 3: Data input sources, Inputs PMSLT and PMSLT mechanisms. The color coding from Data inputs sources to Inputs PMSLT link sources with inputs for the PMSLT. Solid arrows represent final inputs and dashed-arrows represent intermediate inputs that need further processing. Purple coding means a process and green coding represent change in mortality and disability prevalence rates to modify the life table parameters. Black color coding with white color font represent final model outcomes. For both, the life table and disease life table, two sets of each are simulated, one for the baseline and the other for the scenario.

Proportional multi-state life-table model. Three sections are presented in Figure 3: Data input sources, Inputs PMSLT and PMSLT mechanisms. The color coding from Data inputs sources to Inputs PMSLT link sources with inputs for the PMSLT. Solid arrows represent final inputs and dashed-arrows represent intermediate inputs that need further processing. Purple coding means a process and green coding represent change in mortality and disability prevalence rates to modify the life table parameters. Black color coding with white color font represent final model outcomes. For both, the life table and disease life table, two sets of each are simulated, one for the baseline and the other for the scenario.

2.2.2.1 Global Burden of Disease data

First, we explain how to obtain the data, second additional processing to derive data not reported (population) and one-year age groups (original data is in five-year age groups) and last procedure to use Dismod II. Data from the Global Burden of Disease data in Figure 3 can be download from here: http://ghdx.healthdata.org/gbd-results-tool. Figure 4 is a screenshot of the GHDx.

Global Burden of Disease data results tool.

Global Burden of Disease data results tool.

Table 2 specifies the selections to do for each of the tabs in Figure 2.

Once the selections described in Table 2 are made, the option *Download CVS** in the GHDx website is selected. A prompt comes up asking for an email address. The data is sent to the designated email address (within minutes) in ZIP format, unzip and use the code below to read the data (4). Here, we selected data for Greater London and England. The aim is to compare and derive scaling factors as for most cities the data is not available from the GBD and country level data may be used and scaled to the city level. Note that the data input requirement for the PMSLT, except population numbers, is in rates. Therefore the scaling is to better reflect the burden of an area, this is a different issue than working with numbers (e.g. total mortality numbers, total YLDs numbers) as in the ITHIM approach.

  1. Read GBD data
GBDdata <- read.csv("data/UK/gbd2016.csv", stringsAsFactors = F)

The following codes serves to sort out the GBD data to the inputs required for the life table, disease life table and Dismod II

These data should be used to generate the general life table and disease life tables (Figure 3).

  1. Change all upper cases to lower case and delete () from variables.
GBDdata <- mutate_all(GBDdata, funs(tolower))
GBDdata$measure[GBDdata$measure=="ylds (years lived with disability)"] <- "ylds"
  1. Create age categories index in GBDdata (mid age, to match cohort running), total of 20 age groups. These are the age cohorts to simulate.
# 22, 27, 32, 37, 42, 47, 52, 57, 62, 67, 72, 77, 82, 87, 92, 97
GBDdata$age_cat [GBDdata$age =="under 5"] <- 2
GBDdata$age_cat [GBDdata$age =="5 to 9"] <- 7
GBDdata$age_cat [GBDdata$age =="10 to 14"] <- 12
GBDdata$age_cat [GBDdata$age =="15 to 19"] <- 17
GBDdata$age_cat [GBDdata$age =="20 to 24"] <- 22
GBDdata$age_cat [GBDdata$age =="25 to 29"] <- 27
GBDdata$age_cat [GBDdata$age =="30 to 34"] <- 32
GBDdata$age_cat [GBDdata$age =="35 to 39"] <- 37
GBDdata$age_cat [GBDdata$age =="40 to 44"] <- 42
GBDdata$age_cat [GBDdata$age =="45 to 49"] <- 47
GBDdata$age_cat [GBDdata$age =="50 to 54"] <- 52
GBDdata$age_cat [GBDdata$age =="55 to 59"] <- 57
GBDdata$age_cat [GBDdata$age =="60 to 64"] <- 62
GBDdata$age_cat [GBDdata$age =="65 to 69"] <- 67
GBDdata$age_cat [GBDdata$age =="70 to 74"] <- 72
GBDdata$age_cat [GBDdata$age =="75 to 79"] <- 77
GBDdata$age_cat [GBDdata$age =="80 to 84"] <- 82
GBDdata$age_cat [GBDdata$age =="85 to 89"] <- 87
GBDdata$age_cat [GBDdata$age =="90 to 94"] <- 92
GBDdata$age_cat [GBDdata$age =="95 plus"] <- 97
  1. Create age and sex categories to obtain population numbers. Population numbers from GBD are used in Dismod II. For the Life table (Figure 3), the numbers may be from the synthetic population. For now, the Life table is set up with population numbers derived from the GBD data.
GBDdata$sex_age_cat <- paste(GBDdata$sex,GBDdata$age_cat, sep = "_"  )
  1. Convert string variables to numeric to do calculations.
GBDdata$val <- as.numeric(as.character(GBDdata$val))
  1. Generate population numbers for Greater London in a new data frame (“GBD_population”). Note that there is data for England as well, which is used in a separate rmarkdown document (GBDCompare). Population numbers are derived from rates per 100,000 and total numbers of cases.
GBD_population <- filter(GBDdata, measure == "deaths", cause == "all causes", 
metric == "rate" | metric == "number" ) %>% select(metric, age_cat, sex, val, 
sex_age_cat, location)
  1. Generate population numbers from given number of cases and rates per 100,000 people.
 for (i in 1:nrow(GBD_population)) {
    if (GBD_population$metric[i] == "number") {
      GBD_population$five_year_population[i] <- GBD_population$val[i] * 100000/
      GBD_population$val[i + 2]}     
    else {GBD_population$five_year_population[i] <- NA}
 }
  1. Remove rows with zero
GBD_population <- GBD_population[!is.na(GBD_population$five_year_population),]
  1. Keep relevant variables
GBD_population <- filter(GBD_population) %>% 
select(sex_age_cat, sex, age_cat, five_year_population, location)
  1. Create data frames for Greater London to be later used for: a) interpolation of rates and b) PMLT cohorts.
GBD_population_GL <-  filter(GBD_population, location == "greater london") %>% 
select(sex_age_cat, age_cat, sex, five_year_population, location, sex_age_cat)
  1. Check population total numbers
GreaterLondon <- sum(GBD_population_GL$five_year_population)
  1. Generate data frames for Greater London with per person rates (per 100,000 in original data).
GBDGL <- filter(GBDdata, location == "greater london" & metric == "rate") %>% 
select(measure, location, sex, age, metric, cause, val, age_cat, sex_age_cat)
GBDGL$one_rate <- GBDGL$val/100000
  1. For the life table, we need to mortality and total yld rates in one year age intervals. The original data is in five years. Thus, the following code is used to interpolate a single-year age distribution form a five-yearly distribution.
## Loop to generate interpolated rates for all cause mortality and ylds for males and females
## UPDATE WITH YLDS RATES ALL CAUSES ADJUSTED FOR ALL OTHER MODELLED DISEASES.

i_sex <- c("male", "female")
i_measure <- c("deaths", "ylds") #" (years lived with disability)")

for(sex_index in i_sex) {
  for (measure_index in i_measure) {
   
data <- filter(GBDGL, measure == measure_index, sex == sex_index, 
cause == "all causes") %>% select(measure, location, sex, age, metric, 
cause, val, age_cat, one_rate) 
assign(paste(sex_index, measure_index, "interpolated_data", sep = "_"), data)
x=data$age_cat
y=log(data$one_rate)


interpolation_func <- splinefun(x, y, method = "natural", ties = mean)

interpolated <- as.data.frame(interpolation_func(seq(0, 100, 1)))
age <- seq(0,100,by=1)
interpolated <- cbind(interpolated, age)
interpolated[,1] <- exp(interpolated[,1])
colnames(interpolated)[1] <- paste(measure_index)
## Add column with sex to create age_sex category to then merge with input_life table
interpolated$sex <- paste(sex_index)
interpolated$sex_age_cat <- paste(interpolated$sex, interpolated$age, sep = "_")
## Change name of column death to mx and ylds to pyld_rate to then merge 
## with input_life table

if (colnames(interpolated)[1] == "deaths")
    colnames(interpolated)[1] <- paste("mx")
  else
    colnames(interpolated)[1] <- paste("pyld_rate")

# Name data frame
assign(paste(sex_index, measure_index, "interpolated", sep = "_"), interpolated)
# plot(interpolated$age,interpolated$rate)

##Do graph with another layer for original rates in age groups for comparison purposes.
p <- ggplot(data = interpolated,mapping = aes(age, interpolated[,1])) +
              geom_line(aes(color = "Interpolated")) +
              geom_point(
              data = data, 
              mapping = aes(age_cat, one_rate, color = "Original")) +
              labs(colour="",x="Age",y="Rates", sep = " ") +
              labs (title = paste("Rates", sex_index, "all cause", 
              measure_index, sep = " "), size=14) +
              theme_classic() +
              theme (plot.title = element_text(hjust = 0.5))
              print(p)
  }
}
TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

  1. The following code brings together all the required data to create the inputs life table (Figure 3).
## Create empty data frame
input_life_table <- as.data.frame(matrix(0, ncol=2, nrow = 202))
## Give variables names (mx=total mortality rate, 
## yldsx = total years lived with disability rates )
names(input_life_table) <- c("age", "sex")
## Populate life table: Age, 0 to 100
input_life_table[1:101, 1] <- c(0:100)
input_life_table[102:202, 1] <- c(0:100)
## Populate life table: sex, male, female
input_life_table[1:101, 2] <- "male"
input_life_table[102:202, 2] <- "female"
## Create variable age_sex to match with population data. 
input_life_table$sex_age_cat <- paste(input_life_table$sex,input_life_table$age,
                                sep = "_"  )
##Populate life table: five_year_pop. 

input_life_table <- merge(input_life_table, select(GBD_population_GL, c(sex_age_cat, five_year_population)), by = 'sex_age_cat', all = TRUE) 
##Populate life table: mortality rates (from interpolated rates)
## Males mortality rates. NEED MATCHING NAMES
input_life_table <- merge(input_life_table, select(male_deaths_interpolated , c(sex_age_cat, mx)), by = 'sex_age_cat', all.x = TRUE)
input_life_table <- merge(input_life_table, select(female_deaths_interpolated ,c(sex_age_cat, mx)), by = 'sex_age_cat', all.x = TRUE)
input_life_table <- merge(input_life_table, select(male_ylds_interpolated ,c(sex_age_cat, pyld_rate)) , by = 'sex_age_cat', all.x = TRUE) 
input_life_table <- merge(input_life_table, select(female_ylds_interpolated ,c(sex_age_cat, pyld_rate)), by = 'sex_age_cat', all.x = TRUE)

## DISCUSS WITH ALI TO IMPROVE

input_life_table$mx <-  ifelse(input_life_table$sex == "male",
input_life_table$mx.x, input_life_table$mx.y)
input_life_table$pyld_rate <-  ifelse(input_life_table$sex == "male",
input_life_table$pyld_rate.x, input_life_table$pyld_rate.y)
## Drop redundant variables
input_life_table <- subset(input_life_table, select = -c(mx.y, mx.x, pyld_rate.x, pyld_rate.y))

## Cross check population numbers
sum(input_life_table$five_year_population, na.rm = TRUE)
## [1] 8778396
## Sort data by sex and age to use in life table function. 
input_life_table<-input_life_table[order(input_life_table$sex, input_life_table$age),] 
  1. Generate baseline life tables from input_data_frame using run_life_table function. Life tables are generated for each age (5-years) and sex cohort. It is assumed that current observed rates of mortality and disability are going to be observed in the future.
i_age_cohort <- c(22, 27, 32, 37, 42, 47, 52, 57, 62, 67, 72, 77, 82, 87, 92, 97)

general_life_table_list_bl <- list()

index <- 1

for (age in i_age_cohort){
  for (sex in i_sex){
    # cat("age ", age, " and sex ", sex, "\n") #Uncomment to see index
    general_life_table_list_bl[[index]] <- run_life_table(in_idata = input_life_table,
    in_sex = sex, in_mid_age = age)
    index <- index + 1
  }
}

## Uncommnet to check life table list
# View(general_life_table_list_bl[[2]])
  1. Generate baseline disease life tables using run_disase function. Inputs of the disease life tables are: incidence, case fatality (Figure 3, purple shaded area) and disability weights. Internally consistent estimates of incidence and case fatality are derived from GBD data (by five-year age groups and sex) and Dismod II. Disability weights are derived from GBD disease specific YLDs and prevalence (by five-year age groups and sex) adjusted for all-cause YLDs (see section 2.1.1.4).
    Note that Dismod II is an external software. The code below is used to generate the inputs for Dismod II using GBD data (see Figure 3, purple shaded area). Dismod II requires data for the collection (population numbers and all cause mortality rates) and data set (disease specific parameters). We will generate a data set with all inputs for dismod by 5 year-age (represent at mid age) groups and sex for five diseases (ischemic heart disease, ischemic stroke, diabetes mellitus, colon and rectum cancer and breast cancer-women). WE WILL NEED ADJUSTMENT FOR TO MATCH RRS TO DISEASES (COLON AND RECTUM CANCER AND DIABETES MELLITUS). *Code here is work in progress**, for now the disease input data is organised externally with excel.
## Use externaly generated inputs for for disease life table. The previous code will serve to to some of the data preparation for DISMOD. Also need to devlop code to calculated DWs, for now, all generated in excel. 

idata <- read.csv("data/UK/idata.csv", stringsAsFactors = F)
## Use run_disease

i_disease <- c("ihd", "istroke", "diabetes", "colon_cancer", "breast_cancer")
i_sex <- c("male", "female")
i_age_cohort <- c(22, 27, 32, 37, 42, 47, 52, 57, 62, 67, 72, 77, 82, 87, 92, 97)

disease_life_table_list_bl <- list()
index <- 1

for (age in i_age_cohort){
  for (sex in i_sex){
    for (disease in i_disease) {
      # Exclude breast_cancer for Males
      if (sex == "male" && disease == "breast_cancer"){
       # cat("\n") #Uncomment to see list
      }
      else {
        # cat("age ", age, " sex ", sex, "and disease", disease, "\n") #Uncomment to see list
        disease_life_table_list_bl[[index]] <- run_disease(in_idata = idata, in_sex = sex, in_mid_age = age, in_disease = disease)
        index <- index + 1
      }
    }
  }
}
##### Uncommnet to check disease life table list
# View(disease_life_table_list_bl[[8]])
  1. Generate mock change in incidence of disease to generate scenario life tables. In the final version, this will come from the calculated PIFs by disease. PIFs are applied here to incidence, however, may also be applied to case fatality, depending on the RRs. Also, we can include time delays from change in exposure to change in health outcomes via the PIF.
## Create value to use as factor changing incidence rates. 

incidence_change <- 0.95

#####Generate scenario incidence (for each disease)

incidence_sc <- list()
index <- 1

for (age in i_age_cohort){
  for (sex in i_sex){
    for (disease in i_disease) {
      
      # Exclude breast_cancer for Males
      if (sex == "male" && disease == "breast_cancer"){
        # cat("\n") # Uncomment to see list
      }
      else {
        
        incidence_sc[[index]] <- disease_life_table_list_bl[[index]]$incidence_disease * incidence_change
        index <- index + 1
        
      }
    }
  }
}

##### Uncommnet to check scenario incidence
# View(incidence_sc[[1]])
  1. Use scenario incidence to calculate scenario disease life tables.
disease_life_table_list_sc <- list()
index <- 1
for (age in i_age_cohort){
  for (sex in i_sex){
    for (disease in i_disease) {
      # Exclude breast_cancer for Males
      if (sex == "male" && disease == "breast_cancer"){
        # cat("\n")
      }
      else {
        # cat("age ", age, " sex ", sex, "and disease", disease, "\n")
        # modify idata's incidence for the said scenario
        td1 <- idata
        td1[td1$age >= age & td1$sex == sex,][[paste("incidence", disease, sep = "_")]] <- incidence_sc[[index]]
        
        # Instead of idata, feed td to run scenarios
        disease_life_table_list_sc[[index]] <- run_disease(in_idata = td1, in_sex = sex, in_mid_age = age, in_disease = disease)
        disease_life_table_list_sc[[index]]$diff_inc_disease <- disease_life_table_list_sc[[index]]$incidence_disease - disease_life_table_list_bl[[index]]$incidence_disease
        disease_life_table_list_sc[[index]]$diff_prev_disease <- disease_life_table_list_sc[[index]]$px - disease_life_table_list_bl[[index]]$px
        disease_life_table_list_sc[[index]]$diff_mort_disease <- disease_life_table_list_sc[[index]]$mx - disease_life_table_list_bl[[index]]$mx
        disease_life_table_list_sc[[index]]$diff_pylds_disease <- (disease_life_table_list_sc[[index]]$px - disease_life_table_list_bl[[index]]$px) * disease_life_table_list_bl[[index]]$dw_disease
        index <- index + 1
      }
    }
  }
}
##### Uncommnet to check scenario life tables
# View(disease_life_table_list_sc[[1]])
  1. Calculate life tables scenario.
#####Generate total change in mortality rate
##### Sum mortality rate scenarios (mx_sc_total)


mx_sc_total <- list()
l_index <- 1
index <- 1
for (age in i_age_cohort){
  for (sex in i_sex){
    mortality_sum <- NULL
    create_new <- T
    
    for (disease in i_disease) {
      if (sex == "male" && disease == "breast_cancer"){
        # cat("\n")
      }else{
        
        if (create_new){
          mortality_sum <- select(disease_life_table_list_sc[[index]], c('age', 'sex'))
          mortality_sum$total <- 0
          create_new <- F
          mortality_sum$total <- mortality_sum$total + (disease_life_table_list_sc[[index]]$diff_mort_disease)
        }else{
          mortality_sum$total <- mortality_sum$total + (disease_life_table_list_sc[[index]]$diff_mort_disease)
        }
        
        # cat(age, " - ", sex," - ",  disease," - ",  index, " - ", l_index,  "\n")
        index <- index + 1
      }
    }
    mx_sc_total[[l_index]] <- mortality_sum 
    l_index <- l_index + 1
  }
}  

##### Uncommnet to check sceanrio mortality and changes 
# View(mx_sc_total[[1]])


#####Generate total change in prevalent yld rates
#####total ylds rate= sum (change prevalence disease*dw)


pylds_sc_total <- list()
l_index <- 1
index <- 1
for (age in i_age_cohort){
  for (sex in i_sex){
    pylds_sum <- NULL
    create_new <- T
    
    for (disease in i_disease) {
      if (sex == "male" && disease == "breast_cancer"){
        # cat("\n")
      }else{
        
        if (create_new){
          pylds_sum <- select(disease_life_table_list_sc[[index]], c('age', 'sex'))
          pylds_sum$total <- 0
          create_new <- F
          pylds_sum$total <- pylds_sum$total +   (disease_life_table_list_sc[[index]]$diff_pylds_disease)
        }else{
          pylds_sum$total <- pylds_sum$total + (disease_life_table_list_sc[[index]]$diff_pylds_disease)
        }
        
        # cat(age, " - ", sex," - ",  disease," - ",  index, " - ", l_index,  "\n")
        index <- index + 1
      }
      
    }
    pylds_sc_total[[l_index]] <- pylds_sum
    l_index <- l_index + 1
  }
}  

##### Uncommnet to check scenario pyld change
# View(pylds_sc_total[[2]])

##### Calculate general life tables with modified mortality and pylds total#############################
###Original mortality rate is modified by the mx_sc_total (total change in mortality from diseases)
###Original pyld rate is modified by the change in each disease pylds


general_life_table_list_sc <- list()
index <- 1


for (age in i_age_cohort){
  for (sex in i_sex){
    # cat("age ", age, " and sex ", sex, "\n")
    # modify idata's mortality and pyld total for the said scenario
    td2 <- input_life_table
    td2[td2$age >= age & td2$sex == sex,][[paste("mx")]] <- general_life_table_list_bl[[index]]$mx + mx_sc_total[[index]]$total
    td2[td2$age >= age & td2$sex == sex,][[paste("pyld_rate")]] <- general_life_table_list_bl[[index]]$pyld_rate + pylds_sc_total[[index]]$total
    
    # Instead of idata, feed td to run scenarios
    general_life_table_list_sc[[index]] <- run_life_table(in_idata = td2, in_sex = sex, in_mid_age = age)
    
    index <- index + 1
  }
}
##### Uncommnet to check scenario life tables
# View(general_life_table_list_sc[[1]])
# View(general_life_table_list_bl[[2]])
  1. Generate list of outputs by age and sex
#####In the following list "output_life_table", 32 data frames are nested per age and sex cohort

output_burden <- list()
l_index <- 1
index <- 1
for (age in i_age_cohort){
  for (sex in i_sex){
    
    #Males do not have breast cancer, that is why we need the if/else. 
    #We create a TRUE/FALSE variable for the loop to move into the next disease
    
    create_new <- T
    for (disease in i_disease) {
      if (sex == "male" && disease == "breast_cancer"){
        # cat("\n")
      }else{
        
        if (create_new){
          output_burden_sc <- select(disease_life_table_list_sc[[index]], c('age', 'sex', 'incidence_disease', 'mx', 'px'))
          names(output_burden_sc)[names(output_burden_sc) == 'incidence_disease'] <- paste('incidence_disease', disease, "sc", sep = "_")
          names(output_burden_sc)[names(output_burden_sc) == 'mx'] <- paste('mx', disease, "sc", sep = "_")
          names(output_burden_sc)[names(output_burden_sc) == 'px'] <- paste('px', disease, "sc", sep = "_")
          output_burden_bl <- select(disease_life_table_list_bl[[index]], c('incidence_disease', 'mx', 'px'))
          names(output_burden_bl)[names(output_burden_bl) == 'incidence_disease'] <- paste('incidence_disease', disease, "bl", sep = "_")
          names(output_burden_bl)[names(output_burden_bl) == 'mx'] <- paste('mx', disease, "bl", sep = "_")
          names(output_burden_bl)[names(output_burden_bl) == 'px'] <- paste('px', disease, "bl", sep = "_")
          
          ####New list to add calculations for changes in burden of disease (incidence and mortality numbers)
          
          output_burden_change <- list()
          
          output_burden_change$inc_num_bl <- disease_life_table_list_bl[[index]]$incidence_disease * (1 - disease_life_table_list_bl[[index]]$px) * general_life_table_list_bl[[l_index]]$Lx
          output_burden_change$inc_num_sc <- disease_life_table_list_sc[[index]]$incidence_disease * (1 - disease_life_table_list_sc[[index]]$px) * general_life_table_list_sc[[l_index]]$Lx
          output_burden_change$inc_num_diff <- (disease_life_table_list_sc[[index]]$incidence_disease * (1 - disease_life_table_list_sc[[index]]$px) * general_life_table_list_sc[[l_index]]$Lx) - (disease_life_table_list_bl[[index]]$incidence_disease * (1 - disease_life_table_list_bl[[index]]$px) * general_life_table_list_bl[[l_index]]$Lx)
          
          output_burden_change$mx_num_bl <- disease_life_table_list_bl[[index]]$mx * general_life_table_list_bl[[l_index]]$Lx
          output_burden_change$mx_num_sc <- disease_life_table_list_sc[[index]]$mx * general_life_table_list_sc[[l_index]]$Lx
          output_burden_change$mx_num_diff <- (disease_life_table_list_sc[[index]]$mx * general_life_table_list_sc[[l_index]]$Lx) - (disease_life_table_list_bl[[index]]$mx * general_life_table_list_bl[[l_index]]$Lx) 
          
          names(output_burden_change)[names(output_burden_change) == 'inc_num_bl'] <- paste('inc_num_bl', disease, sep = "_")
          names(output_burden_change)[names(output_burden_change) == 'inc_num_sc'] <- paste('inc_num_sc', disease, sep = "_")
          names(output_burden_change)[names(output_burden_change) == 'inc_num_diff'] <- paste('inc_num_diff', disease, sep = "_")
          names(output_burden_change)[names(output_burden_change) == 'mx_num_bl'] <- paste('mx_num_bl', disease, sep = "_")
          names(output_burden_change)[names(output_burden_change) == 'mx_num_sc'] <- paste('mx_num_sc', disease, sep = "_")
          names(output_burden_change)[names(output_burden_change) == 'mx_num_diff'] <- paste('mx_num_diff', disease, sep = "_")
          
          
          ###Bind all lists
          
          output_burden_sc <- cbind(output_burden_sc, output_burden_bl)
          output_burden_sc <- cbind(output_burden_sc, output_burden_change)
          
          create_new <- F
          
          #Here the calculations above are repeated, here is where the F is telling to move into the next disease
          
        }else{
          
          td3 <- select(disease_life_table_list_sc[[index]], c('incidence_disease', 'mx', 'px'))
          names(td3)[names(td3) == 'incidence_disease'] <- paste('incidence_disease', disease, "sc", sep = "_")
          names(td3)[names(td3) == 'mx'] <- paste('mx', disease, "sc", sep = "_")
          names(td3)[names(td3) == 'px'] <- paste('px', disease, "sc", sep = "_")
          
          td4 <- select(disease_life_table_list_bl[[index]], c('incidence_disease', 'mx', 'px'))
          names(td4)[names(td4) == 'incidence_disease'] <- paste('incidence_disease', disease, "bl", sep = "_")
          names(td4)[names(td4) == 'mx'] <- paste('mx', disease, "bl", sep = "_")
          names(td4)[names(td4) == 'px'] <- paste('px', disease, "bl", sep = "_")
          
          
          
          output_burden_change2 <- list()
          
          output_burden_change2$inc_num_bl <- disease_life_table_list_bl[[index]]$incidence_disease * (1 - disease_life_table_list_bl[[index]]$px) * general_life_table_list_bl[[l_index]]$Lx
          output_burden_change2$inc_num_sc <- disease_life_table_list_sc[[index]]$incidence_disease * (1 - disease_life_table_list_sc[[index]]$px) * general_life_table_list_sc[[l_index]]$Lx
          output_burden_change2$inc_num_diff <- (disease_life_table_list_sc[[index]]$incidence_disease * (1 - disease_life_table_list_sc[[index]]$px) * general_life_table_list_sc[[l_index]]$Lx) - (disease_life_table_list_bl[[index]]$incidence_disease * (1 - disease_life_table_list_bl[[index]]$px) * general_life_table_list_bl[[l_index]]$Lx)
          
          output_burden_change2$mx_num_bl <- disease_life_table_list_bl[[index]]$mx * general_life_table_list_bl[[l_index]]$Lx
          output_burden_change2$mx_num_sc <- disease_life_table_list_sc[[index]]$mx * general_life_table_list_sc[[l_index]]$Lx
          output_burden_change2$mx_num_diff <- (disease_life_table_list_sc[[index]]$mx * general_life_table_list_sc[[l_index]]$Lx) - (disease_life_table_list_bl[[index]]$mx * general_life_table_list_bl[[l_index]]$Lx) 
          
          names(output_burden_change2)[names(output_burden_change2) == 'inc_num_bl'] <- paste('inc_num_bl', disease, sep = "_")
          names(output_burden_change2)[names(output_burden_change2) == 'inc_num_sc'] <- paste('inc_num_sc', disease, sep = "_")
          names(output_burden_change2)[names(output_burden_change2) == 'inc_num_diff'] <- paste('inc_num_diff', disease, sep = "_")
          names(output_burden_change2)[names(output_burden_change2) == 'mx_num_bl'] <- paste('mx_num_bl', disease, sep = "_")
          names(output_burden_change2)[names(output_burden_change2) == 'mx_num_sc'] <- paste('mx_num_sc', disease, sep = "_")
          names(output_burden_change2)[names(output_burden_change2) == 'mx_num_diff'] <- paste('mx_num_diff', disease, sep = "_")
          
          
          ###Bind all lists
          
          output_burden_sc <- cbind(output_burden_sc, td3)
          output_burden_sc <- cbind(output_burden_sc, td4)
          output_burden_sc$age_cohort <- age 
          output_burden_sc <- cbind(output_burden_sc, output_burden_change2)
          
        }
        
        # cat(age, " - ", sex," - ",  disease," - ",  index, " - ", l_index,  "\n")
        index <- index + 1
      }
      
    }
    
    # general_life_table_list_sc and general_life_table_list_bl (Lx)
    output_burden_lf_sc <- select(general_life_table_list_sc[[l_index]], c('Lx', 'Lwx'))
    names(output_burden_lf_sc)[names(output_burden_lf_sc) == 'Lx'] <- paste('Lx', "sc", sep = "_")
    names(output_burden_lf_sc)[names(output_burden_lf_sc) == 'Lwx'] <- paste('Lwx', "sc", sep = "_")
    
    output_burden_lf_bl <- select(general_life_table_list_bl[[l_index]], c('Lx', 'Lwx'))
    names(output_burden_lf_bl)[names(output_burden_lf_bl) == 'Lx'] <- paste('Lx', "bl", sep = "_")
    names(output_burden_lf_bl)[names(output_burden_lf_bl) == 'Lwx'] <- paste('Lwx', "bl", sep = "_")
    
    
    output_burden_lf_sc$Lx_diff <- general_life_table_list_bl[[l_index]]$Lx - general_life_table_list_sc[[l_index]]$Lx
    output_burden_lf_sc$Lwx_diff <- general_life_table_list_bl[[l_index]]$Lwx - general_life_table_list_sc[[l_index]]$Lwx
    
    output_burden_sc <- cbind(output_burden_sc, output_burden_lf_sc)
    output_burden_sc <- cbind(output_burden_sc, output_burden_lf_bl)
    
    
    output_burden[[l_index]] <- output_burden_sc
    l_index <- l_index + 1
    
  }
}

#Uncomment to check
# View(output_burden[[32]])
  1. Combine all lists of outputs in a data frame to facilitate extraction of outcomes of interest and plotting.
#####Generate a data frame for all results and create function to get outcomes. 

output_df <- plyr::ldply(output_burden, rbind)

#Remove variables that are not used in the generation of outputs. CHANGE THIS NAMES, TOO LONG

output_df <- subset(output_df, select = -c(incidence_disease_ihd_bl, incidence_disease_ihd_sc, incidence_disease_istroke_bl, incidence_disease_istroke_sc, incidence_disease_diabetes_bl, incidence_disease_diabetes_sc, incidence_disease_breast_cancer_bl, incidence_disease_breast_cancer_sc, 
incidence_disease_colon_cancer_bl, incidence_disease_colon_cancer_sc,  mx_ihd_bl, mx_ihd_sc, mx_istroke_bl, mx_istroke_sc, mx_diabetes_bl, mx_diabetes_sc, mx_breast_cancer_bl, mx_breast_cancer_sc, mx_colon_cancer_bl, mx_colon_cancer_sc, px_ihd_bl, px_ihd_sc, px_istroke_bl, px_istroke_sc, px_diabetes_bl, px_diabetes_sc, px_breast_cancer_bl, px_breast_cancer_sc, px_colon_cancer_bl, px_colon_cancer_sc))
  1. Plot outcomes for each age and sex cohort using created function. ADD A LOOP TO GENERATE ALL OUTCOMES FOR IN_AGE, IN_POPULATION, IN_OUTCOMES AND IN_LEGEND. FIX THE COLORS IN THE FUNCTION.
## Do a loop, the first in_outcomes is age.

i_outcome <- c("mx", "inc")

for (age in i_age_cohort){
  for (sex in i_sex) {
     for (disease in i_disease){
       for (outcome in i_outcome) {
         {
       if (sex == "male" && disease == "breast_cancer"){
         # cat("\n")
       }else{

 p <- plot_output(in_data = output_df, in_age = age, in_population = sex, in_outcomes = c("age", paste(outcome, "num", "bl", disease, sep = "_"), paste(outcome, "num", "sc", disease, sep = "_"), paste(outcome, "num", "diff", disease, sep = "_")), in_legend = paste(outcome, disease, sep = " "))
           }
         }
       }
     }
   }
 }
TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

TRUE

paste(age, sex, disease, outcome, sep = “_“)

  1. Add up all outcomes per year of simulation with generated function.
####Generate data frame with all outputs to graph total change in burden by simulation year. 
###first, need to run function for males and females separetly, in_cohorts indicates the number of age cohorts to 
####include. To show all select 16. what is specified in in_outcomes will be graphed or can also be presented as a total. 
aggregate_frame_males <- gen_aggregate(in_data = output_df, in_cohorts = 16, in_population = "male", in_outcomes = c('inc_num_bl_ihd','inc_num_sc_ihd', "inc_num_diff_ihd"))

aggregate_frame_females <- gen_aggregate(in_data = output_df, in_cohorts = 16, in_population = "female", in_outcomes = c('inc_num_bl_ihd','inc_num_sc_ihd', "inc_num_diff_ihd"))

#####The following adds up both males and females
# Remove non-numeric columns starting with age and sex
aggregate_frame_males <- aggregate_frame_males %>% select(-starts_with("age"), -starts_with("sex"))

aggregate_frame_females <- aggregate_frame_females %>% select(-starts_with("age"), -starts_with("sex"))

# Create a copy of aggregate_frame_females
total_aggr <- aggregate_frame_females
# Add aggregate_frame_males values to it
for (i in 1:ncol(aggregate_frame_females)){
  total_aggr[i] <- total_aggr[i] + aggregate_frame_males[i]
}
total_aggr$sim_year <- seq.int(nrow(total_aggr))
  1. Plots using total aggregated change in burden.
####This plot has to be customised to in_outcomes, here, only totals shown, but specifications are up to the user

####[] is used here to indicate the number of simulation years into the future. 
####Disease outcomes has to be changed to the outcome of interest

disease_outcome <- "Deaths ischemic heart disease"

total_plot <- ggplot(total_aggr[1:20,], aes(x = sim_year)) +
  
  geom_line(mapping = aes(y = total_inc_num_bl_ihd, colour = "total_mx_num_bl_ihd")) +
  theme_classic() +
  geom_hline(yintercept=0, linetype="dashed", color = "black") +
  geom_line(mapping = aes(y = total_inc_num_sc_ihd, colour = "total_mx_num_sc_ihd")) +
  geom_line(mapping = aes(y = total_inc_num_diff_ihd, colour = "total_mx_num_diff_ihd")) +
  xlab ("Simulation years") + ylab ("Cases") + labs (title = paste(disease_outcome)) +
  theme(plot.title = element_text(hjust = 0.5, size = 12)) +
  scale_color_discrete(name = paste(""), labels = c("Baseline", "Difference", "Scenario")) +
  theme(plot.title = element_text(hjust = 0.5))
  

#Print to view
print(total_plot)
TRUE

TRUE

3 Comments

3.1 Road injuries in the PMsLT

The disease model used in each of the disease life table is not directly applicable to road injuries, however, similar concept can be follow. Firstly, changes in road fatalities impact on the overall mortality rate, hence, by knowing the road fatality rates for baseline and scenarios, we will be able to incorporate changes to mortality attributable to road fatalities. For road injuries, methods developed by Kavi Bhalla and Marko Tanio (REFS) that derive the average YLD attributable to life long and short term injuries can be applied to derive the change in total YLDs (CHECK THAT THESE WERE DEVELOPED AS INCIDENCE YLDs).MT’s methods assumes that injuries do not reduce the life expectancy of the injured person.

References

Barendregt, J.J., and J.L. Veerman. 2010. “Categorical Versus Continuous Risk Factors and the Calculation of Potential Impact Fractions.” Journal Article. J Epidemiol Community Health 64 (3): 209–12. doi:10.1136/jech.2009.090274.

Blakely, T., L. J. Cobiac, C. L. Cleghorn, A. L. Pearson, F. S. Deen, G. Kvizhinadze, N. Nghiem, M. McLeod, and N. Wilson. 2015. “Health, Health Inequality, and Cost Impacts of Annual Increases in Tobacco Tax: Multistate Life Table Modeling in New Zealand.” Journal Article. PLoS Med 12. doi:10.1371/journal.pmed.1001856.

Briggs, Adam, Peter Scarborough, and Adrian Smith. 2016. “Modelling in Public Health.” Book Section. In Public Health Intelligence: Issues of Measure and Method, edited by Krishna Regmi and Ivan Gee, 67–90. Cham: Springer International Publishing. doi:10.1007/978-3-319-28326-5_4.

Cobiac, L.J., T. Vos, and J.J. Barendregt. 2009. “Cost-Effectiveness of Interventions to Promote Physical Activity: A Modelling Study.” Journal Article. Plos Med 6 (7): e1000110–e1000110. doi:10.1371/journal.pmed.1000110.

Gold, Marthe R., David Stevenson, and Dennis G. Fryback. 2002. “HALYs and Qalys and Dalys, Oh My: Similarities and Differences in Summary Measures of Population Health.” Journal Article. Annu Rev Public Health 23 (1): 115–34. doi:doi:10.1146/annurev.publhealth.23.100901.140513.

Murray, Christopher J. L., Majid Ezzati, Abraham D. Flaxman, Stephen Lim, Rafael Lozano, Catherine Michaud, Mohsen Naghavi, et al. 2012. “GBD 2010: Design, Definitions, and Metrics.” Journal Article. The Lancet 380 (9859): 2063–6. doi:10.1016/S0140-6736(12)61899-6.

Roux, L., M. Pratt, and T. O. Tengs. 2008. “Cost Effectiveness of Community-Based Physical Activity Interventions.” Journal Article. Am J Prev Med 35. doi:10.1016/j.amepre.2008.06.040.

Vos, T., R. Carter, J. J Barendregt, Mihalopoulos C., JL. Veerman, A. Magnus, L. Cobiac, Bertram MY., and AL. Wallace. 2010. “Assessing Cost-Effectiveness in Prevention (Ace-Prevention): Final Report.” Report.